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We have numerically investigated the dynamics of vortices in a clean layered superconductor 
placed in a perpendicular magnetic field. We describe the energetics using a Ginzburg-Landau free 
energy functional in the lowest Landau level approximation. The dynamics are determined using the 
time-dependent Ginzburg-Landau approximation, and thermal fluctuations are incorporated via a 
Langevin term. The c-axis conductivity at nonzero frequencies, as calculated from the Kubo formal- 
ism, shows a strong but not divergent increase as the melting temperature Tm is approached from 
above, followed by an apparently discontinuous drop at the vortex lattice freezing temperature. The 
discontinuity is consistent with the occurrence of a first-order freezing. The calculated equilibrium 
properties agree with previous Monte Carlo studies using the same Hamiltonian. We briefly discuss 
the possibility of detecting this fluctuation conductivity experimentally. 



I. INTRODUCTION 

Vortices in the mixed state of a clean type-II supercon- 
ductor are believed to break the translational symmetry 
and form a triangular Abrikosov lattice for magnetic field 
B exceeding the lower critical field H c \ . At a sufficiently 
high temperature T, thermal fluctuations are expected 
to melt this lattice and restore the translational sym- 
metry through a solid-liquid phase transition. However, 
in most low-T c superconductors, this melting occurs near 
the upper critical field H C 2, and is thus difficult to distin- 
guish from the usual superconducting-normal transition. 
In the high-T c cuprates, however, this melting transition 
is typically well separated from H c i- 

Many experiments suggest that this solid-liquid phase 
transition is first order. For example, the resistivity 
of untwinned single crystals of YBa2Cu307_5 (YBCO) 
drops sharply at a temperature Tm well below the H C 2 (T) 
linei This temperature also coincides with a discontin- 
uous jump in the magnetization. 2,3 Most distinctively, 
both a latent heat and a specific heat discontinuity have 
been observed at the transition. These signatures have 
been seen in untwinned single crystals of YBCO for mag- 
netic fields both parallel and perpendicular to the c-axis; 4 

The first order nature of the transition is supported 
by a number of theoretical models^SiS^ At least at 
high fields, this transition is thought to represent a si- 
multaneous melting of the vortex lattice in the ab-plane 
and decoupling of the vortex "pancakes" in different lay- 
ers. For example, numerical studies of layered supercon- 
ductors, using Monte Carlo methods applied to a model 
Hamiltonian in the lowest Landau level (LLL) approxi- 
mation, suggest a simultaneous melting and decoupling 
transitioni£i£i! A similar conclusion is also suggested by 
studies of melting using an analogy with a 2D Bose 
system^ On the other hand, some workers have suggested 
that the LLL actually gives no phase transition at all, but 
only a crossover associated with interlayer decoupling. 9 

In this paper, we extend previous numerical studies 
of flux lattice melting to treat the dynamics of a vor- 



tex system. Our calculations are based on a Lawrence- 
Doniach model for the free energy functional of a lay- 
ered superconductor, treated in the LLL approximation, 
and are carried out as a function of temperature at fixed 
magnetic field. The dynamics are treated within the 
time-dependent Ginzburg-Landau approximation, with 
Langevin noise included to simulate the effects of ther- 
mal fluctuations. A previous calculation, for a similar 
model in two dimensions, and with spherical rather than 
periodic boundary conditions, has been carried out by 
Kienappel and Moored. The LLL approximation is ex- 
pected to be most accurate at strong magnetic fields 
(i? C 2/3 < B < H C 2), but may have a slightly broader 
range of validity at low temperatures, since a weak partic- 
ipation of higher Landau levels at such temperatures can 
be incorporated by a suitable renormalization of the LLL 
model parameters^ The LLL approximation fails, how- 
ever, at weak magnetic fields, because it omits the effects 
of thermally induced vortex-antivortex pairs. Consistent 
with previous Monte Carlo studies, we find a single first- 
order liquid-solid phase transition with simultaneous loss 
of in-plane and interplane vortex order. However, the 
Langevin simulation also yields information about dy- 
namical properties such as the conductivity. In particu- 
lar, we find that the c-axis conductivity shows a striking, 
but not divergent, increase as the first-order melting tem- 
perature Tm{B) is approached from above. 

The remainder of this paper is organized as follows. In 
Sec. II we describe the Langevin model, and our method 
for calculating various static and dynamic quantities from 
the model. Our results are presented in Sec. Ill, followed 
by a brief discussion in Sec. IV. 



II. FORMALISM 

A. Model Hamiltonian and Dynamical Equations 

We consider a three-dimensional (3D) superconductor 
consisting of a stack of Josephson- or proximity-coupled 
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2D layers. We assume that this system is described by 
the free energy functional 



(1) 



Here n is the layer index, d is the thickness of one layer, 
and 

/»hMr)] = a(T)|^„(r)| 2 + i/3|^(r)| 4 
-iKV-^) V„(r) 
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ipn(r) is the order parameter of the nth layer, g = — 2e is 
the charge of a Cooper pair, e? is the distance between the 
layers, and ct(T), (3 , m a b, and m c are material-dependent 

parameters. The phase factor Xn,n+i = /id +1 ' )<i dzA z , 
where $o = hc/2e is the flux quantum, and A is the 
vector potential. We will assume that the external mag- 
netic field B|jz, i.e. perpendicular to the layers, so that 
Xn,n+i = 0, and we choose a gauge such that A = —Byx. 
We also neglect screening currents and fluctuations of the 
vector potential, so that the local and externally applied 
magnetic fields are equal. This should be a good approx- 
imation when the Ginzburg-Landau parameter k > 1, as 
in the cuprate superconductors. Finally, we assume that 
B is uniform throughout the superconductor. 

In the LLL approximation, the order parameter in each 
layer is expanded as 



■0„(r) = ip ^2 c n,k<t>k(x, y), 
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where 
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are the lowest eigenstates of the kinetic energy operator 
(— iTiV — q A/c) 2 /(2m* b ), corresponding to eigenvalue 

hqB/(2m ab c). Here ^ = [iri 2 \a H (T)| 2 /(£ 2 (3 2 )] 1/4 , 
a H (T) = a(T)(l - B/H c2 ), I = {\q\B /Tic)- 1 / 2 is the 
magnetic length, and t Q = (4tt/ v / 3) 1/2 L The magnitude 
of ipo is chosen so that the spatial average of | V^tt. ( r ) 1 2 1S 
\<Xh\/{P@a) with /3a = 1.169... when the vortices are 
arranged in a triangular lattice. 

We assume the system is a parallelepiped of dimen- 
sions L x , L y , and L z , with periodic boundary conditions 
in all three directions. We choose L z — N z d, where 
N z is an integer. The periodicity condition in the x- 
direction ip n (x + L x ,y) = ip n (x,y) implies k = 2i:m/L x , 
where m is an integer. If each layer contains vor- 
tices, there will be N$ independent c n .k's labeled by 
m = 0, . . . , (N^ — 1). The periodicity constraint in 
the y-direction, \ip(x, y + L y )\ — \ip(x,y)\, implies that 
Cn m' = c n m f° r au m ' — m modulo Nj,. Finally, 



the periodicity constraint in the z-direction implies that 
C n+N z ,k — Cn,k- Besides these periodicity conditions, the 
cell dimensions in the x- and y-direction are chosen to be 
compatible with a possible triangular lattice. This choice 
may be written as L x /L y — 2n x /(v3n 3/ ) where n x and 
n y are the number of vortices along a given row or column 
parallel to the x- or y-direction, and = n x n y . 

Using Eq. © for the order parameter, we can rewrite 
Eq. as: 



r( n ) I t(") 
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where 1S ^ ne ^ ree ener §y P er layer, and is the 
coupling between the nth and (n + l)th layers. These 
terms take the form 



^d/W) = g 2 (B,T)n x \sgn(a H )Y, 



\Cn,k\ 



^n,k ^n,k-\-p ^n,k-\-q ^n,fc+p 



k,p,q 



and 



F { c ] /{k B T) = g 2 (B, T) n xV J2 |c*,« - c k , n+1 \ 2 . (7) 
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Here, we have defined 



;(p, q) = y&/^exp [-1 2 (p 2 + q 2 ) /2] 



g 2 (B,T) = irfdoa 2 H /(f3k B T), 
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and introduced the dimensionless interlayer coupling 
strength n = J/\an\ where J = H 2 /(2m c d 2 ) is the 
Josephson coupling between the layers. The quantity 
sgn(a^f) = —1 or +1 in the superconducting or normal 
regimes; the mean field instability occurs when an(T) = 
0. 

The parameter g 2 (B,T) represents the ratio of the 
superconducting condensation energy per vortex per 
layer -Kl 2 doa 2 H / '/3 to the thermal energy k B T within the 
Ginzburg-Landau approximation. Note that, for fixed 
a B , g 2 (B,T) varies inversely with temperature. Thus, 
a plot of system properties as a function of |<?|, may be 
viewed as a plot as a function of T for fixed B; however, 
small \g\ represents large T (vortex liquid phase). Previ- 
ous calculations^ have provided evidence that there is a 
first-order vortex lattice melting transition as a function 
oig 2 . 

We study the dynamics of this system using the time 
dependent Ginzburg-Landau (TDGL) equation in the 
presence of a Langevin noise term. We write this equa- 
tion as 



, dip n (T, t) 
dt 
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+ fn(r,«), 



(10) 
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where F is the relaxation time parameter, and £(r,i) is 
a white noise term characterized by the correlation func- 
tions 

(£»(r,t)) = 
(CMKn'(r',0)) = 2k B TTS(r-r')5 n ^S(t), 

where (• • •) denotes an ensemble average. We assume 
that r is real because the system has particle-hole sym- 
metry; with our choice of units, T has the same dimen- 
sions as h. The noise term ensures that the system will 
remain in a steady state at temperature T. 

Langevin dynamical calculations have previously been 
carried out by Ryu and Stroud 12 to study vortex lat- 
tice melting for both clean and dirty high-T c layered su- 
perconductors. They differ from the present calculation 
by using a different equilibrium free energy functional T 
than ours. In the model of Ref. Il2t the flux lines can- 
not be broken; this feature should lead to rather different 
results from those obtained in the present model calcu- 
lations. 

In the present LLL expansion, the TDGL equation may 
be rewritten as 



dc n 



dr 



sgn(a H ) Cn,k 



+ 2 H v ^ k ~ p > q ^ c « 
- V (Cn-l,k - 2c„ 



n,p-\-q n,k-\-p n,q 



'71+1 



k)+£'n,fc(l-)(ll) 



Here we have introduced a dimensionless time variable 
t = |a#|£/r = tj to, where to = T/\au\ is a character- 
istic relaxation time. The noise term £' is now described 
by the correlation functions^ 



(e n Ar)) = 0; 



(£'n,/c( T ) C'n'.fe'( T ')> 



(12) 



B. Calculated Quantities 

1. Equilibrium Quantities 

Equilibrium quantities can be computed as time aver- 
ages of the solutions to the TDGL equations, either in 
the solid or the liquid phase. According to the ergodic 
hypothesis, this procedure should give the same results 
as an equilibrium average obtained by treating Eq. (jSJ) 
as a Hamiltonian. Wc have, in fact, confirmed this point 
by comparing some of our results with those obtained 
earlier by other workers from Eq. JSJ using Monte Carlo 
techniques £1 

We have evaluated several thermodynamic properties 
of the system. One is the generalized Abrikosov factor 



(3 A = N z L x L y 



(14) 



When the vortices form a triangular lattice, /3a reaches 
its minimum value of /3a , but exceeds this value for other 
vortex configurations. We have also computed the spatial 
average of \t/j n (r)\ 2 , defined as 
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Tab 



N z L x L y 



]T J d 2 i#„(r)| 2 , 



(15) 



which at low temperatures reaches the mean field value 
r ofc F = \ a H \/(0aP)- Both Pa and r a b vary smoothly with 
temperature and thus do not show any special behavior 
at the flux lattice melting temperature Tm(B). We have 
therefore also examined three other equilibrium quanti- 
ties which more clearly show signals of this transition: 
the isothermal shear modulus n(T) of the flux lattice; a 
quantity we denote C(T), which measures the degree of 
coherence between vortices in adjacent layers; and the zz- 
component of the helicity modulus tensor, T CC (T), which 
measures stiffness against a long-wavelength twist in the 
phase of the order parameter. 

The shear modulus n(T) is defined^ by 
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(d 2 F 



N z L x L y V d6 2 



(16) 
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where 6 is the shear angle. The free energy F can be 
obtained from Eq. JSJ), using F = —ksTlnZ where 
Z = Tr e _JC 'l ksT and the trace is taken over the classical 
variables Cn,k and c* k . An explicit 2D form for a(T) in 
the LLL approximation has been given in Ref. Il4l where 
it has been found that ^t(T) reaches its mean field value 
M MF (T) = 0.354A^fc B 7y(T) at low T and it vanishes 
everywhere in the liquid phase. If the transition between 
the vortex solid and vortex liquid state is first-order, then 
/x(T) will, in the thermodynamic limit, jump discontin- 
uously from a finite value to zero at Tm(B). However, 
such a jump does not prove that the melting transition 
is first-order, since certain continuous melting transitions 
in two dimensions also have a jump in /i(T) at melting. 15 
However, other independent calculations give evidence 
that the melting transition is first-order within the LLL 
approximation in three dimensions (e.g., by exhibiting a 
finite latent heat). 
C(T) is defined by 



C(T) 



E»/rf 2 r|y; n+1 (r)-^(r)| ; 
2E„/rf 2 r|^(r)| 2 



(17) 



At low T, where the vortex system forms a flux lattice 
with flux lines all parallel to the c-axis, ipn+i(r) — V'n(r) 
and hence C = 0. By contrast, deep in the liquid phase, 
the phases of ^> n (r) and ip n+ \(r) are uncorrected, and C 
approaches unity. To calculate C(T) and other ratios of 
spatial averages, we evaluate the ratios at fixed time, and 
then average over a period of time as described below. 

Finally, the helicity modulus component T CC (T) is de- 
fined by the relation^ 



T, 



m = - 



1 fd 2 F 



da 2 , 



(18) 
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Here a z is an additional uniform vector potential ap- 
plied in the z-direction (besides that which is needed to 
produce the magnetic field), and V is the system vol- 
ume. A further discussion of the meaning of T cc is to 
be found in Ref. In the mean field approximation, 
T CC (T) is approximated by T™ F (T) = 2 Jdd r^ /<^ 2 , 
where r^ F = \cxh\/(PaP) is the mean field value of the 
quantity r a b defined in Eq. (|15|) . T cc is shown in Ref. 
to drop discontinuously to zero at Tm(B). 



2. Dynamical Quantities 

The wave-number- and frequency-dependent conduc- 
tivity of the vortex system can be computed using the 
Kubo formula. If the frequencies lu satisfy the condition 
hu> -C ksT, one may use the Kubo formula in the classi- 
cal limit >ii 



o>/(q,w) 



1 



dt I d 3 xd 3 x' 



k B TV , 

elq .(x-x')--* (jAi(Xji)>(x ' ;0)) . (19) 



Here (x, t) is the /xth component of the current den- 
sity, and (T^i/ (q, u>) is the [iv th component of the complex 
conductivity tensor for a wave number q and frequency 
ui, and (• • •) denotes an average over the thermal noise 
distribution. 

In the present work, we have considered only the con- 
ductivity component o~ c , which requires only the c-axis 
current density. Within the Lawrence-Doniach model, 
this current density, for z in the region between the nth 
and (n + l)th layer, is 



h q 

m r d 



lm[r n+1 (r)Mr)}- (20) 



If we expand ip(r) using the representation we find 
that J z (t) = w: J2 n Jd 2 rji n) (r,t) is 



m c d 



where 



k°n,k ■ 



(21) 



(22) 



Note that this current density includes only the Joseph- 
son currents between the layers, and not any additional 
normal currents which may be flowing in parallel. 

The corresponding real fluctuation conductivity 
a c .i(ui) = Re[cr C c(q = 0,tt>)] follows from the Kubo for- 
mula lfT§j) : 



d 2 N 2 



k B TV 



dtcos(ojt)(J z (t)J z (0)). (23) 



Upon using Eq. (|21|l . it takes the form 
a cA (cj') N z n x G 2 (T) r°° 



era 



b y JO 



dTCOs(wV) (J(r)J(O)), 



where cto = q tolV'ol / m c — w ^Oj and 
2 ,_, 2doV9 2 (B,T) 



G 2 (T) 



(24) 



(25) 



Besides the frequency-dependent conductivity, it is 
sometimes of interest to compute the integrated fluctu- 
ation conductivity, 72, defined byi£ 



1 f°° 

72 = - / dw a- c i(u). 
n Jo 



With the use of Eq. (|24|) , 72 can be simplified to 



72 



G 2 (T)a n x 
n y N z 



(\J(o)\ 2 ) 



(26) 



(27) 



where J(0) is given by Eq. (|2*2*|l . 



III. RESULTS 

We have solved the Langevin equations (jl 1|> numeri- 
cally using a second order Runge-Kutta algorithm. In 
this algorithm, T is correct through C(e 2 ) in the time 
step e»i*i In most of our simulations, we used e = 0.15 to; 
a smaller time step of e = 0.05 to was found to give simi- 
lar results but to require more computer time. The real 
and imaginary parts of the noise term £' n k (r) in Eq. l(TT|) 
are chosen from Gaussian distributions with a mean zero 
and a variance cr 2 /e where a 1 = 2/[n x g 2 (B, T)]. This 
choice insures that these terms have mean and variance 
which satisfy Eqs. (12) and (13). 

In most cases, we have started our simulations from 
the low-temperature Abrikosov phase, then gradually in- 
creased T, taking the initial state for a higher T as the 
equilibrium state for the previous slightly lower T. We 
have verified that our results exhibit only a little hys- 
teresis - that is, we obtain the same equilibrium and 
nearly the same dynamical results, whether Tm(B) is 
approached from below or from above. We have found 
that our choice of initial state generally has little effect on 
dynamics, providing we "anneal" our sample for a long 
enough time as described in the next paragraph. We have 
confirmed this lack of effect by obtaining similar results 
for various calculated dynamical quantities whether we 
begin by choosing an Abrikosov or a liquid-like initial 
state. 

For each temperature considered, we have allowed the 
system to equilibrate for a period ranging from 10 3 to 
4 x 10 6 time steps, before starting to compute averages, 
the larger number corresponding to temperatures close 
to Tm- We then run the dynamics for an additional 3 x 
10 4 to 10 6 time steps at this temperature, and use these 
results to compute the averages. 
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FIG. 1: (a) Ratio of the mean-square gap r a b [Eq. 1151 1 to 
its mean-field value r^ F , and the ratio of the generalized 
Abrikosov factor /3a [Eq. I|14jl ] to its value /3a in a triangular 
lattice, plotted as a function of \g\ [Eq. ©]. (b) The calcu- 
lated ratios of the shear modulus /i and the c-axis helicity 
modulus T cc to their mean- field values fiMF and Y,?f F , as 
obtained from Langevin dynamical simulations, plotted as a 
function of | g \ . Also shown is C(T) [Eq. 1171 1. The lattice used 
contains N$ = 6x6 vortices in the ab-plane and N z — 12 lay- 
ers in the c-direction, with periodic boundary conditions; the 
interlayer coupling is chosen so that r]\g\ = 0.02. The error 
bars in this and later figures represent the standard deviations 
of results from about five Langevin dynamical simulations run 
for equal lengths of time. 



To calculate the quantity of interest, we include in the 
averages only the results obtained in every Nq time steps, 
where Nq is chosen as explained below. If this procedure 
is used, then, according to Ref.|2(J, the consecutive values 
included in the average become nearly statistically inde- 
pendent. We choose Nq using a criterion involving the 
so-called self-correlator. This self-correlator is defined by 
the relation C x (k) = ((x i+k Xi) - (x l ) 2 )/((xf) - (xi) 2 ), 
where Xi is the physical quantity of interest at the ith. 
time step, and (...) is a time average. With this defini- 
tion, C X (Q) = 1; also, C x (k) decreases as k increases. The 
optimum choice of k to be used in the simulations (i.e., 
the optimum number of steps between those included in 



the averages) is that which makes C x (k) as small as possi- 
ble, typically less than 0.05. For our simulations, we find 
that this optimum value of k = iVo is typically between 
20 and 50. In general, we find that collective properties 
such as conductivity require much longer runs than sin- 
gle particle properties such as r ab \ the optimum ratio of 
collective to single-particle running times itself depends 
on the system size. 

As a test of our numerical algorithm, we have com- 
puted several equilibrium properties of the system which 
have been previously evaluated using Monte Carlo 
methods, 6,7 Our results are shown in Fig. ^ For these 
calculations, the lattice used contained 6x6 vortices in 
the ab-plane and 12 layers in the c-direction. The cou- 
pling between the layers is chosen as rj\g\ — 0.02. We 
emphasize that our calculations are intended to probe a 
range of physically reasonable parameters, rather than to 
describe any specific superconductor. In Ref. 0, f]\g\ was 
estimated to vary from 0.0075 (BSCCO) to 0.30 (YBCO) 
in typical cuprate superconductors at temperatures and 
magnetic fields where the LLL approximation is likely to 
be valid; our choice falls well within this range. 

Figure nja) displays the generalized Abrikosov factor 
Pa/ Pa [Eq. O] and the quantity r ab [Eq. (I5)l] as a 
function of g [Eq. @]. These results are similar to previ- 
ous results obtained using the Monte Carlo methodi 7 - In 
Fig. IHb) we show the calculated /V/x MF [Eq. JSJ] ver- 
sus \g\. The sharp drop in n/n MF near \g\ » 4 is clearly 
visible. Although our sample sizes are quite small, the 
calculated equilibrium quantities still show the expected 
behavior of larger samples, though somewhat broadened 
by the substantial finite size effects. 

Also shown in Fig. ^b) is the interlayer coupling 
strength parameter C(T) [Eq. I|17|) ] and the helicity mod- 
ulus component T CC (T), both plotted as functions of \g\. 
Both of these quantities (which are sensitive to z-axis co- 
herence) show a drop near \g\ w 4, but C(T) is expected 
to vary smoothly through this region of \g\, while T CC (T) 
is expected to drop discontinuously to zero in the ther- 
modynamic limit; some evidence of this distinction can 
be seen in the Figure. The simultaneous drop in fj,//i MF 
and T CC (T) near Tm suggests that there is simultaneous 
flux lattice melting in the ab-plane and interlayer decou- 
pling in the c-direction, near T = Tm, consistent with 
previous calculations in clean systems 

Next, we turn to the dynamics of the system. We have 
calculated (J z (t/ to)J z (0)} as a function of various system 
parameters. In Fig. [21 we plot this correlation function as 
a function of tj to for several values of \g\ both above and 
below the expected melting point, denoted |<7jw|. (|<?m| ~ 
4 for lattices of this size and our choice of parameters^* 7 -) . 
Clearly, the decay rate slows considerably as melting is 
approached from higher temperatures, i.e., from smaller 
values of \g\. However, the decay rate is rapid and only 
weakly dependent on |<?| in the vortex lattice phase, \g\ > 
\9m\- 

To make this |g|-dependence more apparent, we plot in 
Fig-Elthe half-life Ti/ 2 of this correlation function versus 
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FIG. 2: Normalized correlation function Cj(r) = {J'(t)J'(0)}/{J'(0) 2 ) plotted as a function of r for several values of g as 
indicated in the Figure. The freezing transition occurs at \g\ « 4 for this lattice size. Note the expanded horizontal scale. The 
lattice used here has 16 vortices in each plane, and 16 planes (16 x 16 lattice), and we choose r]\g\ = 0.05. 



\g\ for two different system sizes, normalized by N z to. 
Ti/2 is denned as the time at which (J z (t/ to)J z (0)) has 
fallen to half its r = value. Consistent with Fig. 
7]/2 is relatively small in both the solid and the liquid 
phases far from |<?m|; it increases noticeably as \gtd\ is 
approached from the liquid but not from the solid phase. 
Despite this increase, we believe that T\ /2 will not diverge 
at \qm\i because \gu\ corresponds to a first-order melting 
transition, with no dynamical critical phenomena such as 
a diverging correlation time. 

Figures0fa) and EJb) show cr C) i(a/), as obtained from 
Eq. 1)23(1 for several values of ui', and for two different 
system sizes. The lines simply connect the calculated 
points. The melting value \gn\ ~ 4, as estimated from 
Fig. 1(b). For uj' = 0.004 and 0.008, a c ,i(u') increases 
strongly as |(/m| is approached from the liquid side. There 
is a smaller increase at higher a/, because the transition 
primarily affects fluctuations on a lower frequency scale. 
In the solid phase, there is little evidence of fluctuations 
in <7 C) i(u/), which remains very small at nonzero frequen- 
cies for all \g\ > \gu\ studied. At fixed |gr| in the liquid 
phase near \gn\, we expect <7 C) i(u/) to decrease mono- 
tonically with increasing u)'\ we ascribe any deviation 
from monotonic behavior in Fig. 4 is to numerical un- 
certainties. Similarly, although <7 c ,i(a/) sometimes seems 
to peak at a value of \g\ slightly smaller than \gu\i we 
believe that this behavior also lies within our numerical 
uncertainties. 

In Fig. [3] we show the total integrated fluctuation con- 
ductivity 72 [Eq. l|77jl ]. in units of a G 2 (T), plotted as a 
function of \g\ for two lattice sizes. We have computed 72 
using the equilibrium expression Eq. (|27|l , which is equiv- 
alent to the frequency integral of the quantity shown in 
Fig. 0Ja) or Bib). As \ghi\ is approached from the liq- 




\g\ 

FIG. 3: Half-life T1/2 characterizing the decay rate of 
{Jz{t/ to)Jz(O)) plotted as a function of \g\. The full lines 
simply connect the calculated points. Error bars have the 
same meaning as in Fig. The sizes of the lattice used are 
indicated in the Figure legend and the interlayer coupling is 
chosen so that rj\g\ — 0.05. 



uid side; 72 falls sharply at |<7m|, to a small value in 
the vortex lattice phase. This drops is expected to be 
discontinuous in the thermodynamic (large size) limit. 
As previously, the full lines simply connect the calcu- 
lated points. Note that 72/ [cr G 2 (T)] is approximately 
|g|-independent in the liquid phase, because of the way 
it is normalized (G 2 oc g 2 ). 
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FIG. 4: (a) The real part of the fluctuation conductivity, 
<7 c ,i(a/)/[<7oG 2 (T)], plotted versus Igl for several values of uj 
for a 4 x 4 x 10 = 16 x 10 lattice, (b) Same as (a) except that 
the lattice size is 4x4x 16. Note the sharp increase in er Cil (o/) 
in the vortex liquid phase near and below \givi\ for small a/. 
The parameters are the same as in Fig. |2] For clarity, we do 
not show error bars; they are comparable to those of Fig. 3. 



IV. DISCUSSION 

The present results are consistent with the scenario of a 
first-order melting transition at Tm(B) 7 where long-range 
vortex order in the ab-planes, and in the c-direction, dis- 
appear simultaneously. This interpretation is supported 
by the behavior of the hclicity modulus T CC (T) and shear 
modulus (J-(T), both of which vanish at the same tem- 
perature T. Although the present calculations are lim- 
ited to relatively small samples (with fewer than 50 vor- 
tex pancakes per plane), similar behavior has been ob- 
served in Monte Carlo simulations for considerably larger 
systems 

The behavior of dynamical properties, such as er c 
is also consistent with first-order melting. For small u>', 
cr C) i(a/) shows a strong increase as T approaches Tm from 
above. This behavior occurs because, in the solid phase, 
the vortex pancakes in adjacent layers lie above one an- 
other; as a result, fluctuating currents between the layers 



FIG. 5: The integrated fluctuation conductivity 

72/[o"oG 2 (T)] plotted as a function of \g\ for two lattice 
sizes as shown in the Figure legend. The parameters are the 
same as in Fig. [5] 



are small. On the other hand, when the lattice melts (low 
Ifl 1 ! or high T), the vortex pancakes in adjacent layers no 
longer lie directly above one another; hence, fluctuating 
phase gradients between the layers increase the current 
fluctuations and the fluctuation conductivity. In the liq- 
uid phase, er c i(ij/) decreases with decreasing \g\, because 
Tj/2 is becoming smaller. However, cr cl (w') once again 
appears to show a discontinuity rather than a divergence 
at | gA/ 1, consistent with the first-order nature of the tran- 
sition. 

Why does cr c i(a/) in the liquid state, shown in Fig. 0] 
for several frequencies, decrease with increasing temper- 
ature T above Tm? We believe this decrease occurs be- 
cause, as shown in Fig. [31 T\j2 decreases with increasing 
T. By contrast, the quantity j2/(&oG 2 ), shown in Fig.0 
behaves differently from a Ci \\ it has a discontinuity at Tm 
but varies slowly in the liquid. The lack of any clear peak 
in 72/(coG 2 ) near Tm can be understood from Eq. 126|l . 
which shows that 72 is independent of T1/2, depending 
only on equal-time current density fluctuations at t = 0. 
By contrast, cr c i(a/) is sensitive to T1/2 especially for 
small u>'. 

At this point, we briefly comment on a seemingly coun- 
terintuitive feature of the dynamical results shown in 
Figs.0J namely, that a Cj i(oj') is small in the vortex lattice 
phase for nonzero a/, even well below Tm ■ Intuitively, one 
might expect, since this phase is superconducting, with 
a finite helicity modulus in the c direction, that a c .\{uj') 
would be large in this regime. However, this behavior is 
actually physically reasonable; our picture of the under- 
lying physics is the following. We believe that a Ct ±(u>') 
corresponding to our model dynamics is the sum of two 
parts: (i) the fluctuation conductivity shown in Figs.0Ja) 
and Efb) (whose integral is shown in Fig. |5j and (ii) a 
delta function at zero frequency, corresponding to per- 
fect conductivity. The delta function does not appear in 
Fig. Ufa) or Qfb) because those calculations are carried 



8 



out at finite frequency, nor does it appear in the integral 
shown in Fig. [5] The strength of this delta function is 
proportional to the helicity modulus shown in Fig.^b), 
which vanishes for T > Tm- Thus, o~ Ct x(u)') is small be- 
low Tm simply because the fluctuation contributions to 
the conductivity are small in this temperature range; the 
system is still phase coherent in the c-direction and still 
has a finite helicity modulus below Tm- Although it may 
seem strange that a c .i(u>') is small for T < Tm and finite 
u/, this behavior is not unprecedented. For example, in 
low-T c s-wave superconductors, the existence of a finite 
gap below T c means that o~i(cu) = for T < T c and for 
huj smaller than twice the energy gap. 

To our knowledge, no direct measurements of a c ,i{u') 
have been carried out in the cuprate superconductors in 
the high-field, clean-limit regime where our calculations 
might be applicable. We therefore comment briefly on 
an entirely different experiment in which the reported 
behavior somewhat resembles that shown in Fig.0] This 
is a recent study of the frequency-dependent conductivity 
of BiSr 2 Ca 2 Cu08 +( 5 within the ab plane at zero applied 
magnetic field>2i This experiment reports a rather sharp 
peak near T c in the real part of the in-plane conductivity 
at about 0.2 THz. This peak is thought to be due to 
fluctuations in the phase of the order parameter which 
are strongest near T c , and weaker both above and below 
T c . We believe that similar fluctuations (probably in the 
amplitude of the order parameter as well as the phase) are 
producing the increase in o~ Ct i{u>') i n the present model 
near Tm- These fluctuations are, we believe, limited in 
size because the melting transition is first-order rather 
than continuous, and they are relatively small for T < 
Tm- 

The present calculations may be relevant to c-axis 
transport at strong magnetic fields (where the LLL ap- 
proximation is adequate) in a clean high-T c supercon- 
ductor (where a first-order vortex lattice melting is ex- 
pected), provided a Langevin dynamics is appropriate. 
We have tried to estimate the numerical value of our cal- 
culated <7 Ci i(u/) [Eq. JUJ] for reasonable experimental 
parameters. For co' <C 1 or uto <C 1, <7 c ,i(u/) has the 
same order of magnitude as the quantity 

a G 2 (T) = q 2 t (2d /%£ 2 n 2 )k B T(rig) 2 g 2 . (28) 

All the quantities in this equation are easily determined 
except to. We attempt to estimate to using an early pa- 
per by Schmid,2£ in which the time-dependent Ginzburg- 
Landau equation is derived from the original BCS the- 
ory within the Gor'kov approximation. Schmid finds 
that t m V[32fc B T c0 (l - T/T c0 )(l - B/H c2 )}, where 
T c o is the mean-field superconducting transition temper- 
ature at B = 0. Taking T c0 ~ 80 K, T ~ 60 K we 



find t ~ 4 x 10" 14 sec at a field of 5 T. Substi- 
tuting this value into Eq. (|28|) . and using ijg = 0.05, 
g 2 ~ 20, and d ~ 5 A, we obtain ooG 2 ~ 5 x 10 11 esu, or 
about 0.06 f2 _1 cm -1 . This conductivity is considerably 
smaller than the apparent c-axis conductivity in the vor- 
tex liquid state, even in a very anisotropic material such 
as BiSr2Ca2CuOg+5 ^ Thus, it might be difficult to ob- 
serve the fluctuation contribution in a clean anisotropic 
superconductor, unless we have substantially underesti- 
mated to- Such an underestimate is possible, since the 
calculations of Schmid are based on a microscopic the- 
ory which may not be directly applicable to the layered 
high-T c materials. 

Few experiments appear to have measured the c-axis 
resistivity at the high fields where the LLL approximation 
would be most accurate. Fuchs et al 23 working at far 
lower fields (~ 25-200 Oe), have observed a simultaneous 
disappearance of resistivity in the ab-plane and in the c- 
direction (indicating a single phase transition), and an 
abrupt increase in the c-axis resistivity at a temperature 
just above that transition. However, their experiments 
are done at such frequencies (~ 72 Hz) that the inductive 
contribution is very dominant in the solid phase. 

Finally, we briefly discuss the fact that our calculated 
hysteretic effects are very small in the vicinity of the first- 
order melting transition. In a real experiment, one might 
expect some evidence of superheating or supercooling. 
This minimal amount of hysteresis may be due to the 
long annealing time before we begin to calculate thermo- 
dynamic averages. Because of this long annealing, our 
system can apparently attain its thermodynamic state of 
minimum free energy before we start computing averages. 

In summary, we have studied both the equilibrium and 
the dynamical behavior of a layered superconductor in 
a strong magnetic fields by solving the time-dependent 
Ginzburg-Landau equations, in the lowest Landau level 
approximation. The effects of fluctuations are incorpo- 
rated by means of a Langevin noise term. The equilib- 
rium properties are found to exhibit behavior similar that 
found in previous Monte Carlo results &§iLii a first-order 
melting transition of the vortex lattice, with a simultane- 
ous loss of in-plane and interplane order. The dynamical 
properties show a strong, but not divergent, increase in 
the c-axis conductivity as Tm is approached from above, 
with a corresponding increase, but no divergence, in the 
half-life T]y 2 of the c-axis current fluctuations. 
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